## Copyright (C) 1996, 1997 John W. Eaton
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2, or (at your option)
## any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, write to the Free
## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
## 02110-1301, USA.

## -*- texinfo -*-
## @deftypefn {Function File} {} fftfilt (@var{b}, @var{x}, @var{n})
##
## With two arguments, @code{fftfilt} filters @var{x} with the FIR filter
## @var{b} using the FFT.
##
## Given the optional third argument, @var{n}, @code{fftfilt} uses the
## overlap-add method to filter @var{x} with @var{b} using an N-point FFT.
##
## If @var{x} is a matrix, filter each column of the matrix.
## @} # deftypefn

## Author: Kurt Hornik <Kurt.Hornik@ci.tuwien.ac.at>
## Created: 3 September 1994
## Adapted-By: jwe

   
FftFilter <- function(b, n) {
  res = list(b = b, n = n)
  class(res) = "FftFilter"
  res
}
 
Filter.FftFilter <- function(filt, x, ...) 
  fftfilt(filt$b, x, filt$n)
  
fftfilt  <- function(b, x, n = NULL) {
  fftfiltx  <- function(b, x, n = NULL) {
 
    ## If n is not specified explicitly, we do not use the overlap-add
    ## method at all because loops are really slow.  Otherwise, we only
    ## ensure that the number of points in the FFT is the smallest power
    ## of two larger than N and length(b).  This could result in length
    ## one blocks, but if the user knows better ...
  
    ## TO DO: fix to work with multicolumn x, so don't have to use apply
    ## below.
    
    N = n
    
    l_x = NROW(x)
    l_b = length(b)
   
    if (is.null(n)) {
      ## Use FFT with the smallest power of 2 which is >= length (x) +
      ## length (b) - 1 as number of points ...
      N = 2 ^ (ceiling(log(l_x + l_b - 1) / log(2)))
      B = Fft(postpad(b, N))
      y = ifft(Fft(postpad(x, N)) * B)
    } else {
      ## Use overlap-add method ...
      if (length(n) > 1)
        stop("fftfilt: n has to be a scalar")
      N = 2 ^ (ceiling(log(max(N, l_b)) / log(2)))
      L = N - l_b + 1
      B = Fft(postpad(b, N))
      #B = B[,rep(1., l_x)]
      R = ceiling(l_x / L)
      y = array(0, l_x)
      for (r  in  1:R) {
        lo = (r - 1) * L + 1
        hi = min(r * L, l_x)
        tmp = array(0, N)
        tmp[1:(hi-lo+1)] = x[lo:hi]
        tmp = ifft(Fft(postpad(tmp, N)) * B)
        hi  = min(lo+N-1, l_x)
        y[lo:hi] = y[lo:hi] + tmp[1:(hi-lo+1)]
      }
    }
    y = y[1:l_x]
   
    ## Final cleanups: if both x and b are real respectively integer, y
    ## should also be
   
    if (is.numeric(b) && is.numeric(x)) 
      y = Re(y)
    if (!any(b - round(b))) {
      idx = !any(x - round(x))
      y[idx] = round(y[idx])
    } 
    y
  }

  if (NCOL(x) > 1) {
    apply(x, function(x) fftfiltx(b, x, n = n))
  } else {
    fftfiltx(b, x, n = n)
  }
}

